/**
 * @file   routine_linear_response.c
 * @author xlhu <xlhu@Nano>
 * @date   Thu May 23 19:55:21 2013
 * 
 * @brief  
 * Subroutines for linear response of CDDFT
 * 
 */

#include "routine_linear_response.h"

/**
 * 1. load the ground state of molecule from data file
 * 
 */

void load_ground_state(char* file_name, int cutstate, int st_num,//where? Howmany?
		       double *e_values, double *occ_num,  // orbital informations
		       double *rho_0, double *wfs,         // density and orbitals
		       double *v_0_ext, double *v_0_H, double *v_0_xc,
		       int GsNx, int GsNy, int GsNz,
		       int Nx, int Ny, int Nz) 
{
  int ix=0, iy=0, iz=0;
  int sx=0, sy=0, sz=0;
  char file[80];
  FILE *f_in;

  int state=0;  
  double temp_discard, *temp;
  double offset=1.0e-3;

  sx = (GsNx-1)/(Nx-1);
  sy = (GsNy-1)/(Ny-1);
  sz = (GsNz-1)/(Nz-1);
  printf("*** sx = %d, sy =%d, sz =%d ***\n", sx,sy,sz);
  
  /********
	   1. wfs
  *********/ 
  temp = (double *)malloc(GsNx*GsNy*GsNz*sizeof(double));
  sprintf(file, "%s/wfs.txt\0", file_name);
  f_in = fopen(file, "r");
  //ignore first #cutstate states
  for(state=0;state<cutstate;state++){  
    for(iy=0;iy<GsNy;iy++){
      for(ix=0;ix<GsNx;ix++){
	for(iz=0;iz<GsNz;iz++){  
	  fscanf(f_in,"%lf",&temp_discard); 
	}
      }
    }
  }
  //keep the rest states 
  for(state = 0; state < st_num; state++){ 
    // read into temporal memory
    for(iy=0;iy<GsNy;iy++){
      for(ix=0;ix<GsNx;ix++){
	for(iz=0;iz<GsNz;iz++){
	  fscanf(f_in,"%lf", &temp[iz+GsNz*(ix+GsNx*iy)]); 
	}
      }
    }
    // and pick out the necessary one
    for(iy=0;iy<Ny;iy++){
      for(ix=0;ix<Nx;ix++){
	for(iz=0;iz<Nz;iz++){
	  wfs[iz+Nz*(ix+Nx*(iy+Ny*state))] = temp[iz*sz+GsNz*(ix*sx+GsNx*iy*sy)];
	}
      }
    }
  } 
  fclose(f_in);
  
  /*******
	  2. rho_0
  ********/
  sprintf(file, "%s/rho_0.txt\0", file_name);
  f_in = fopen(file, "r");
  for(iy=0;iy<GsNy;iy++){
    for(ix=0;ix<GsNx;ix++){
      for(iz=0;iz<GsNz;iz++){
	fscanf(f_in,"%lf",&temp[iz+GsNz*(ix+GsNx*iy)]); 
	temp[iz+GsNz*(ix+GsNx*iy)]=fabs(temp[iz+GsNz*(ix+GsNx*iy)])+offset;
      }
    }
  }  
  for(iy=0;iy<Ny;iy++){
    for(ix=0;ix<Nx;ix++){
      for(iz=0;iz<Nz;iz++){
	rho_0[iz+Nz*(ix+Nx*iy)]=temp[iz*sz+GsNz*(ix*sx+GsNx*iy*sy)];
      }
    }
  }      
  fclose(f_in);

  /*******
	  3. v_0_ext
  ********/  
  sprintf(file, "%s/v_ext.txt\0", file_name);
  f_in = fopen(file, "r");
  for(iy=0;iy<GsNy;iy++){
    for(ix=0;ix<GsNx;ix++){
      for(iz=0;iz<GsNz;iz++){
	//fscanf(fout1,"%lf",&v_0_ext[iz+Nz*(ix+Nx*iy)]);  
	fscanf(f_in,"%lf",&temp[iz+GsNz*(ix+GsNx*iy)]); 
      }
    }
  }    
  for(iy=0;iy<Ny;iy++){
    for(ix=0;ix<Nx;ix++){
      for(iz=0;iz<Nz;iz++){
	v_0_ext[iz+Nz*(ix+Nx*iy)]=temp[iz*sz+GsNz*(ix*sx+GsNx*iy*sy)];
      }
    }
  }       
  fclose(f_in); 

  /*******
	  4. v_0_H
  ********/
  sprintf(file, "%s/v_H.txt\0", file_name);
  f_in = fopen(file, "r");
  for(iy=0;iy<GsNy;iy++){
    for(ix=0;ix<GsNx;ix++){
      for(iz=0;iz<GsNz;iz++){
	//fscanf(fout1,"%lf",&v_0_H[iz+Nz*(ix+Nx*iy)]);  
	fscanf(f_in,"%lf",&temp[iz+GsNz*(ix+GsNx*iy)]);  
      }
    }
  }      
  for(iy=0;iy<Ny;iy++){
    for(ix=0;ix<Nx;ix++){
      for(iz=0;iz<Nz;iz++){
	v_0_H[iz+Nz*(ix+Nx*iy)]=temp[iz*sz+GsNz*(ix*sx+GsNx*iy*sy)];
      }
    }
  }      
  fclose(f_in);  

  /*******
	  5. v_0_xc
  ********/
  sprintf(file, "%s/v_xc.txt\0", file_name);
  f_in = fopen(file, "r");
  for(iy=0;iy<GsNy;iy++){
    for(ix=0;ix<GsNx;ix++){
      for(iz=0;iz<GsNz;iz++){
	//fscanf(fout1,"%lf",&v_0_xc[iz+Nz*(ix+Nx*iy)]);  
	fscanf(f_in,"%lf",&temp[iz+GsNz*(ix+GsNx*iy)]);   
      }
    }
  }      
  for(iy=0;iy<Ny;iy++){
    for(ix=0;ix<Nx;ix++){
      for(iz=0;iz<Nz;iz++){
	v_0_xc[iz+Nz*(ix+Nx*iy)]=temp[iz*sz+GsNz*(ix*sx+GsNx*iy*sy)];
      }
    }
  }             
  fclose(f_in);   

  /*******
	  6. occ_num
  ********/
  sprintf(file, "%s/occ_num.txt\0", file_name);
  f_in = fopen(file, "r");
  for(ix=0;ix<cutstate;ix++){
    fscanf(f_in,"%lf",&temp_discard); 
  }
  for(ix=0;ix<st_num;ix++){
    fscanf(f_in,"%lf",&occ_num[ix]); 
  } 
  fclose(f_in); 
  
  /*******
	  7. e_values
  ********/
  sprintf(file, "%s/e_values.txt\0", file_name);
  f_in = fopen(file, "r");
  for(ix=0;ix<cutstate;ix++){
    fscanf(f_in,"%lf",&temp_discard); 
  }
  for(ix=0;ix<st_num;ix++){
    fscanf(f_in,"%lf",&e_values[ix]); 
  } 
  fclose(f_in);    

  // free the temporal memory to read functions
  free(temp);

}


/**
 *        2. evaluate rho and J in subdomain 
 * 
 */
void eval_rho_j(double complex *del_rho, double complex *del_j_1, 
		double complex *del_j_2, double complex *del_j_3, // return values
		double complex *P, double complex omega, 
		double *wfs, double *gradient_wfs, 
		double *e_values, double *occ_num, 
		int N_sample, int occ, int unocc, int Nx, int Ny, int Nz)
{
  int ix,iy,iz,isample, pos;
  int iocc, iunocc, pos_j, pos_b;

  double complex n_np_r,  j_inter;
  double complex n_dnp_1, n_dnp_2, n_dnp_3;
  double complex np_dn_1, np_dn_2, np_dn_3;

  double complex *rho_inter, *j1_inter, *j2_inter, *j3_inter;
  rho_inter = (double complex *)malloc(N_sample*sizeof(double complex));
  j1_inter = (double complex *)malloc(N_sample*sizeof(double complex));
  j2_inter = (double complex *)malloc(N_sample*sizeof(double complex));
  j3_inter = (double complex *)malloc(N_sample*sizeof(double complex));

  double complex p_val;
  // rho and J
  for(ix=0;ix<Nx;ix++){
    for(iy=0;iy<Ny;iy++){
      for(iz=0;iz<Nz;iz++){
	for(isample=0;isample<N_sample;isample++){
	  rho_inter[isample]=0; 
	  j1_inter[isample]=0;    
	  j2_inter[isample]=0;   
	  j3_inter[isample]=0;               
	}
	for(iocc=0;iocc<occ;iocc++){
	  pos_j = iz+Nz*(ix+Nx*(iy+Ny*iocc));
	  for(iunocc=0;iunocc<unocc;iunocc++){
	    pos_b = iz+Nz*(ix+Nx*(iy+Ny*(iunocc+occ)));

	    n_np_r = wfs[pos_j]*wfs[pos_b];  

	    n_dnp_1 = wfs[pos_j]*gradient_wfs[0 + 3*pos_b];
	    n_dnp_2 = wfs[pos_j]*gradient_wfs[1 + 3*pos_b];
	    n_dnp_3 = wfs[pos_j]*gradient_wfs[2 + 3*pos_b];

	    np_dn_1 = wfs[pos_b]*gradient_wfs[0 + 3*pos_j];
	    np_dn_2 = wfs[pos_b]*gradient_wfs[1 + 3*pos_j];
	    np_dn_3 = wfs[pos_b]*gradient_wfs[2 + 3*pos_j];

	    for(isample = 0; isample < N_sample; isample++){
	      p_val = P[iunocc + unocc*(iocc + occ*isample)];

	      rho_inter[isample] += -p_val * n_np_r*occ_num[iocc];
	      j_inter = p_val*omega*occ_num[iocc]/     // two "-" makes "+"
		(e_values[iocc] - e_values[iunocc+occ]); 
 
	      j1_inter[isample] += (-csqrt(-1.0)/2.0)*(n_dnp_1 - np_dn_1)*j_inter;
	      j2_inter[isample] += (-csqrt(-1.0)/2.0)*(n_dnp_2 - np_dn_2)*j_inter; 
	      j3_inter[isample] += (-csqrt(-1.0)/2.0)*(n_dnp_3 - np_dn_3)*j_inter; 
	    } 
	  }
	} 
	for(isample = 0; isample < N_sample; isample++){
	  pos = iz + Nz*(ix + Nx*(iy + Ny*isample));
	  del_rho[pos] = rho_inter[isample]; 
	  del_j_1[pos] = j1_inter[isample]; 
	  del_j_2[pos] = j2_inter[isample];
	  del_j_3[pos] = j3_inter[isample];
	}

      }  //iz
    }  //ix
  }  // iy

  free(rho_inter); free(j1_inter); free(j2_inter); free(j3_inter);

}


        
/**
 * 2. calculate the system matrix, 
 * which mainly based on various potentials: 
 *     del_v_ks, del_A_EM, del_A_xc
 * for orbital pair(iocc, iunocc) of sample l 
 * via (3.2) and (3.3)
 *
 */
void update_matrix(double complex *S0, double complex *S1,
		   double complex *AEM_1,   double complex *AEM_2, 
		   double complex *AEM_3,           // EM field
		   double complex *del_rho, double complex *del_j_1, 
		   double complex *del_j_2, double complex *del_j_3,   
		   double complex omega,  double *rho_0,
		   double *drho_0_x,double *drho_0_y,double *drho_0_z,
		   double *wfs, double *grad_wfs,
		   double *e_values, double *occ_num,  
		   double complex *fxc_kernel, double complex *fxcL0, 
		   double complex *fxcL, double complex *fxcT,
		   double complex *eta_xc, double complex *xi_xc,
		   double complex *gradient_eta_xc, 
		   double complex *hessian_eta_xc,     // help functionals 
		   int N_sample, int occ, int unocc,  //2.Quantum config 
		   int Nx, int Ny, int Nz,         // 3.computational mesh
		   double hx, double hy, double hz,
		   double *Intl_coef,
		   double *QD_emass, double *QD_epsilon)
{
  int i,j,k,l, ix, iy, iz, idx,pos,pos_j, pos_b;
  int iocc, iunocc, isample;
  int st_num = occ + unocc;  // the used orbitals in total
 
  int Msize = N_sample*occ*unocc;
  int Nsize = Nx*Ny*Nz;
  
  // local variables
  double complex np_n_r;
  double complex np_dn_1, np_dn_2, np_dn_3;
  double complex n_dnp_1, n_dnp_2, n_dnp_3;
  double complex dn_1_dnp_2, dn_1_dnp_3, dn_2_dnp_3;
  double complex dnp_1_dn_2, dnp_1_dn_3, dnp_2_dn_3;

  double complex u_1, u_2, u_3, div_u, curl_u_1, curl_u_2, curl_u_3;
  double complex v_1, v_2, v_3, div_v, curl_v_1, curl_v_2, curl_v_3;

  double complex del_v, I1, I2, I3, I4, I5, I6;
  double complex u1, u2, u3, u4, u5, u6, u7, u8;

  double complex *temp_del_rho;
  double complex *temp_del_j_1, *temp_del_j_2, *temp_del_j_3;

  // 1. integration of potentials
  double complex Iv0_ks, IA0_em, IA0_xc;
  double complex *temp_Iv_ks, *temp_IA_EM, *temp_IA_xc; // space 
  double complex *Intl_Iv_ks, *Intl_IA_EM, *Intl_IA_xc; // orbital
  temp_Iv_ks = (double complex *)malloc(Nsize*sizeof(double complex));
  temp_IA_EM = (double complex *)malloc(Nsize*sizeof(double complex));
  temp_IA_xc = (double complex *)malloc(Nsize*sizeof(double complex));
  Intl_Iv_ks = (double complex *)malloc(Msize*sizeof(double complex));
  Intl_IA_EM = (double complex *)malloc(Msize*sizeof(double complex));
  Intl_IA_xc = (double complex *)malloc(Msize*sizeof(double complex));

  // 2. jacobian of del_j
  double complex *dj_1_x, *dj_1_y, *dj_1_z;
  dj_1_x = (double complex *)malloc(Nsize*sizeof(double complex));
  dj_1_y = (double complex *)malloc(Nsize*sizeof(double complex));
  dj_1_z = (double complex *)malloc(Nsize*sizeof(double complex));
  double complex *dj_2_x, *dj_2_y, *dj_2_z;
  dj_2_x = (double complex *)malloc(Nsize*sizeof(double complex));
  dj_2_y = (double complex *)malloc(Nsize*sizeof(double complex));
  dj_2_z = (double complex *)malloc(Nsize*sizeof(double complex));
  double complex *dj_3_x, *dj_3_y, *dj_3_z;
  dj_3_x = (double complex *)malloc(Nsize*sizeof(double complex));
  dj_3_y = (double complex *)malloc(Nsize*sizeof(double complex));
  dj_3_z = (double complex *)malloc(Nsize*sizeof(double complex));

  // 3. Hatree poitential
  double complex *v_H, *rho;   
  v_H = (double complex *)malloc(Nsize*sizeof(double complex));
  rho=(double complex *)malloc(Nsize*sizeof(double complex));

  /**
   * step 1. run over each sample molecule  
   * to evaluate the integration of various potentials 
   *
   */
  for(isample = 0; isample < N_sample; isample++){

    /* pick out the del_rho and del_j for a single sample */
    temp_del_rho = del_rho + Nsize*isample;
    temp_del_j_1 = del_j_1 + Nsize*isample;
    temp_del_j_2 = del_j_2 + Nsize*isample;
    temp_del_j_3 = del_j_3 + Nsize*isample;
    
    /* Call poisson solver for v_H */  
    for(i = 0; i < Nsize; i++) rho[i] = temp_del_rho[i]/QD_epsilon[i];
    poissonFFT(v_H, rho, Nx, Ny, Nz, hx, hy, hz);

    /* Get the gradient of del_j */          
    gradient_of_complex_function(temp_del_j_1, dj_1_x, dj_1_y, dj_1_z, 
				 Nx, Ny, Nz, hx, hy, hz);
    gradient_of_complex_function(temp_del_j_2, dj_2_x, dj_2_y, dj_2_z, 
				 Nx, Ny, Nz, hx, hy, hz);
    gradient_of_complex_function(temp_del_j_3, dj_3_x, dj_3_y, dj_3_z, 
				 Nx, Ny, Nz, hx, hy, hz);

    /* evaluate the poitential for all positions */
    for(iocc = 0; iocc < occ; iocc++){
      for(iunocc = 0; iunocc < unocc; iunocc++){  

	for(ix = 0; ix < Nx; ix++){
	  for(iy = 0; iy < Ny; iy++){
	    for(iz = 0; iz < Nz; iz++){ 

	      // 1. 	      
	      pos_j = iz+Nz*(ix+Nx*(iy+Ny*iocc));// pos of occ orbital
	      pos_b = iz+Nz*(ix+Nx*(iy+Ny*(iunocc+occ))); // pos of unocc orb
	      pos = iz + Nz*(ix+Nx*iy);   // 
		  
	      // 2. p follow b
	      np_n_r = wfs[pos_b]*wfs[pos_j];
	      // d follow grad
	      n_dnp_1 = wfs[pos_j]*grad_wfs[0 + 3*pos_b];
	      n_dnp_2 = wfs[pos_j]*grad_wfs[1 + 3*pos_b];
	      n_dnp_3 = wfs[pos_j]*grad_wfs[2 + 3*pos_b];  	
	      np_dn_1 = wfs[pos_b]*grad_wfs[0 + 3*pos_j];
	      np_dn_2 = wfs[pos_b]*grad_wfs[1 + 3*pos_j];
	      np_dn_3 = wfs[pos_b]*grad_wfs[2 + 3*pos_j];
	      //
	      dn_1_dnp_2 = grad_wfs[0 + 3*pos_j]*grad_wfs[1 + 3*pos_b];
	      dn_1_dnp_3 = grad_wfs[0 + 3*pos_j]*grad_wfs[2 + 3*pos_b];
	      dn_2_dnp_3 = grad_wfs[1 + 3*pos_j]*grad_wfs[2 + 3*pos_b];
	      dnp_1_dn_2 = grad_wfs[0 + 3*pos_b]*grad_wfs[1 + 3*pos_j];
	      dnp_1_dn_3 = grad_wfs[0 + 3*pos_b]*grad_wfs[2 + 3*pos_j];
	      dnp_2_dn_3 = grad_wfs[1 + 3*pos_b]*grad_wfs[2 + 3*pos_j];

	      // 3.1 value of v_ks
	      del_v = v_H[pos] + fxc_kernel[pos]*temp_del_rho[pos]; 
	      temp_Iv_ks[pos] = np_n_r*del_v;

	      // 3.2 value of A_EM
	      temp_IA_EM[pos] = ((np_dn_1 - n_dnp_1)*AEM_1[pos + Nsize*isample] + 
				 (np_dn_2 - n_dnp_2)*AEM_2[pos + Nsize*isample] + 
				 (np_dn_3 - n_dnp_3)*AEM_3[pos + Nsize*isample])*
		(-csqrt(-1.0)/2.0)/c_vacc/QD_emass[0];

	      // 3.3 value of A_xc
	      // u
	      u_1 = temp_del_j_1[pos]/rho_0[pos];
	      u_2 = temp_del_j_2[pos]/rho_0[pos];
	      u_3 = temp_del_j_3[pos]/rho_0[pos];
	      // 
	      div_u = (csqrt(-1.0)*omega*temp_del_rho[pos] - 
		       (drho_0_x[pos]*u_1 + 
			drho_0_y[pos]*u_2 + 
			drho_0_z[pos]*u_3))/rho_0[pos]; 
	      // 
	      curl_u_1 = ((dj_3_y[pos] - dj_2_z[pos]) + 
			  (u_2*drho_0_z[pos] - u_3*drho_0_y[pos]))/rho_0[pos];
	      curl_u_2 = ((dj_1_z[pos] - dj_3_x[pos]) + 
			  (u_3*drho_0_x[pos] - u_1*drho_0_z[pos]))/rho_0[pos];
	      curl_u_3 = ((dj_2_x[pos] - dj_1_y[pos]) + 
			  (u_1*drho_0_y[pos] - u_2*drho_0_x[pos]))/rho_0[pos];
	      // v	      
	      v_1 = -csqrt(-1.0)/2.0*(np_dn_1 - n_dnp_1)/rho_0[pos];
	      v_2 = -csqrt(-1.0)/2.0*(np_dn_2 - n_dnp_2)/rho_0[pos];
	      v_3 = -csqrt(-1.0)/2.0*(np_dn_3 - n_dnp_3)/rho_0[pos];
	      // 
	      div_v = (csqrt(-1.0)*(e_values[iunocc+occ] - e_values[iocc])*
		       wfs[pos_j]*wfs[pos_b] - 
		       (drho_0_x[pos]*v_1 + 
			drho_0_y[pos]*v_2 + 
			drho_0_z[pos]*v_3))/rho_0[pos];
	      // 
	      curl_v_1 = (- csqrt(-1.0)*(dn_2_dnp_3 - dnp_2_dn_3) + 
			  (v_2*drho_0_z[pos] - v_3*drho_0_y[pos]))/rho_0[pos];
	      curl_v_2 = (- csqrt(-1.0)*(dnp_1_dn_3 - dn_1_dnp_3) + 
			  (v_3*drho_0_x[pos] - v_1*drho_0_z[pos]))/rho_0[pos];
	      curl_v_3 = (- csqrt(-1.0)*(dn_1_dnp_2 - dnp_1_dn_2) + 
			  (v_1*drho_0_y[pos] - v_2*drho_0_x[pos]))/rho_0[pos];
	      // 
	      I1 = eta_xc[pos]*(curl_v_1*curl_u_1 + 
				curl_v_2*curl_u_2 + 
				curl_v_3*curl_u_3);
	      I2 = div_v*eta_xc[pos]*div_u;
	      I3 = div_v*(gradient_eta_xc[0 + 3*pos]*u_1 + 
			  gradient_eta_xc[1 + 3*pos]*u_2 + 
			  gradient_eta_xc[2 + 3*pos]*u_3);
	      I4 = (v_1*(hessian_eta_xc[0 + 3*(0 + 3*pos)]*u_1 + 
			 hessian_eta_xc[1 + 3*(0 + 3*pos)]*u_2 + 
			 hessian_eta_xc[2 + 3*(0 + 3*pos)]*u_3)+
		    v_2*(hessian_eta_xc[0 + 3*(1 + 3*pos)]*u_1 + 
			 hessian_eta_xc[1 + 3*(1 + 3*pos)]*u_2 + 
			 hessian_eta_xc[2 + 3*(1 + 3*pos)]*u_3)+
		    v_3*(hessian_eta_xc[0 + 3*(2 + 3*pos)]*u_1 + 
			 hessian_eta_xc[1 + 3*(2 + 3*pos)]*u_2 + 
			 hessian_eta_xc[2 + 3*(2 + 3*pos)]*u_3));
	      I5 = (v_1*gradient_eta_xc[0 + 3*pos] + 
		    v_2*gradient_eta_xc[1 + 3*pos] + 
		    v_3*gradient_eta_xc[2 + 3*pos])*div_u;
	      I6 = div_v*xi_xc[pos]*div_u;

	      // Finally
	      temp_IA_xc[pos] = (I1 + 4.0/3.0*I2 + 2.0*I3 + 2.0*I4 + 2.0*I5 + I6)*
		                csqrt(-1.0)/omega/QD_epsilon[pos]/QD_emass[0];
 
	    }  // iz
	  }  // iy
	}  // ix

	/* Do the numverical integration*/
	Iv0_ks = 0; IA0_em = 0; IA0_xc = 0; 
	for(k = 0; k < Nsize; k++){
	  Iv0_ks += Intl_coef[k]*temp_Iv_ks[k];
	  IA0_em += Intl_coef[k]*temp_IA_EM[k];
	  IA0_xc += Intl_coef[k]*temp_IA_xc[k];
	}

	idx = iunocc + unocc*(iocc + occ*isample);  //position
	Intl_Iv_ks[idx] = Iv0_ks*hx*hy*hz*0.125;    
	Intl_IA_EM[idx] = IA0_em*hx*hy*hz*0.125;
	Intl_IA_xc[idx] = IA0_xc*hx*hy*hz*0.125;

      }  // iunocc
    }  // iocc
  } // isample

  /**
   * step 2. setup basic vectors for executing S*p
   * 
   */
  for(isample = 0; isample < N_sample; isample++){
    for(iocc = 0; iocc < occ; iocc++){
      for(iunocc = 0; iunocc < unocc; iunocc++){
	idx = iunocc + unocc*(iocc + occ*isample);
	S0[idx] = (e_values[iocc] - e_values[iunocc + occ])*
	  Intl_Iv_ks[idx] + 2*omega*(Intl_IA_EM[idx] + Intl_IA_xc[idx]);
	S1[idx] = (e_values[iocc] - e_values[iunocc + occ])*
	  (e_values[iocc] - e_values[iunocc + occ]) - omega*omega;
      }
    }
  }

  /**
   * Finally, free the memory spaces
   * 
   */
  free(temp_IA_EM); free(temp_IA_xc); free(temp_Iv_ks); free(rho);
  free(Intl_IA_xc); free(Intl_IA_EM); free(Intl_Iv_ks); free(v_H); 

  free(dj_1_x);free(dj_1_y);free(dj_1_z);
  free(dj_2_x);free(dj_2_y);free(dj_2_z);
  free(dj_3_x);free(dj_3_y);free(dj_3_z);    
}


/**
 *  update the right hand side of Casida's equation
 *  based on the induced and incident light
 */

void update_rhs(double complex *F,             // output
		double complex *AEM_1, 
		double complex *AEM_2, 
		double complex *AEM_3, 
		double complex omega,
		double *wfs, double *grad_wfs, 
		int N_sample, int occ, int unocc, 
		int Nx, int Ny, int Nz,
		double hx, double hy, double hz,
		double *Intl_coef, double *QD_emass)
{
  int i, k, ix, iy, iz, idx, pos_j, pos_b;
  int iocc, iunocc, isample;

  double complex np_dn_1, np_dn_2, np_dn_3;
  double complex n_dnp_1, n_dnp_2, n_dnp_3;

  int Msize = N_sample*occ*unocc;
  int Nsize = Ny*Nx*Nz;
  double complex int_val, *IA, *Intl_IA_0;
  IA = (double complex *) malloc(Nsize*N_sample*sizeof(double complex));
  Intl_IA_0  = (double complex *)malloc(Msize*sizeof(double complex));

  /// Compute Intl_A_0
  for(iocc = 0; iocc < occ; iocc++){ 
    for(iunocc = 0; iunocc < unocc; iunocc++){
      // 1.
      for(isample = 0; isample < N_sample; isample++) IA[isample] = 0;
      // 2.
      for(iy=0;iy<Ny;iy++){
	for(ix=0;ix<Nx;ix++){
	  for(iz=0;iz<Nz;iz++){   
	    pos_j = iz+Nz*(ix+Nx*(iy+Ny*iocc)); 
	    pos_b = iz+Nz*(ix+Nx*(iy+Ny*(iunocc+occ)));
	    
	    n_dnp_1 = wfs[pos_j]*grad_wfs[0 + 3*pos_b];
	    n_dnp_2 = wfs[pos_j]*grad_wfs[1 + 3*pos_b];
	    n_dnp_3 = wfs[pos_j]*grad_wfs[2 + 3*pos_b];        
	    
	    np_dn_1 = wfs[pos_b]*grad_wfs[0 + 3*pos_j];
	    np_dn_2 = wfs[pos_b]*grad_wfs[1 + 3*pos_j];
	    np_dn_3 = wfs[pos_b]*grad_wfs[2 + 3*pos_j];
	    
	    for(isample = 0; isample < N_sample; isample++){
	      idx = iz + Nz*(ix + Nx*(iy + Ny*isample));
	      IA[idx] = ((n_dnp_1 - np_dn_1)*AEM_1[idx] +
			 (n_dnp_2 - np_dn_2)*AEM_2[idx] +
			 (n_dnp_3 - np_dn_3)*AEM_3[idx])*
		(-csqrt(-1.0))/2.0/c_vacc;
	    }
	    
	  }  //iz
	}  //ix
      } //iy
      // 3. 
      for(isample=0;isample<N_sample;isample++){
	idx = iunocc + unocc*(iocc + occ*isample);
	int_val = 0.0;
	for(k = 0; k < Nsize; k++) int_val += Intl_coef[k]*IA[isample*Nsize+k];
	Intl_IA_0[idx] = int_val*0.125*hx*hy*hz/QD_emass[0];
      } 
      
    } // iocc
  } // iunocc
  
  // then get the right hand side
  for(i = 0; i < Msize; i++) F[i] = -2.0*omega/c_vacc*Intl_IA_0[i];

  //
  free(Intl_IA_0); free(IA);
}


/**
 * 7. BiCGSTAB subroutines for solving Casida's equation
 *    based on TD-CDFT, without update rho and j and AEM during iteration
 * almost the same as for other sparse matrix,
 * The main differnce is that the matrix S is not explicitly known
 * However, S*x = S1.*x + S0
 */

int BiCGStab_Casida(double complex *S0, double complex *S1,
		    double complex *x, double complex *b,
		    int Msize, int precond, int MAXIT, double TOL)
{
  int i,j,k,l, nverbose = MAXIT/10;

  double normb, normr, resid;
  double complex temp, temp1, temp2;

  double complex rho_0, rho_1;
  double complex alpha, beta, w;

  double complex *r, *r0, *p, *phat, *s, *shat, *t, *v;
  r    = (double complex *) malloc(Msize*sizeof(double complex));
  r0   = (double complex *) malloc(Msize*sizeof(double complex));
  p    = (double complex *) malloc(Msize*sizeof(double complex));
  phat = (double complex *) malloc(Msize*sizeof(double complex));
  s    = (double complex *) malloc(Msize*sizeof(double complex));
  shat = (double complex *) malloc(Msize*sizeof(double complex));
  t    = (double complex *) malloc(Msize*sizeof(double complex));
  v    = (double complex *) malloc(Msize*sizeof(double complex));

  /** step 1: r0 = b - Ax0 = b **/
  for(temp1 = 0.0, temp2 = 0.0, k = 0; k < Msize; k++){
    temp1 += conj(b[k])*b[k];
    
    //    temp = (Sp2[k] - k2)*x[k] + Sp1[k]*Intl_v[k] + 2*omega*Intl_A[k]; 
    r[k] = b[k] - (S1[k]*x[k] + S0[k]);
    temp2 += conj(r[k])*r[k];
  }
  normb = sqrt(cabs(temp1)); 
  normr = sqrt(cabs(temp2));
  resid = normr/normb;
  printf("\t ||r_0|| = %20.10lf \n", normr);
  if (resid <= TOL) {
    printf("Good initial guess of x.\n");
    free(r); free(r0); free(p); free(phat); 
    free(s); free(shat); free(t); free(v);
    return 0;
  }

  /** step2: choose r0 fixed **/      
  for(k = 0; k < Msize; k++) r0[k] = r[k];

  /** step3: **/
  rho_0 = 1.0; alpha = 1.0; w = 1.0;

  /** step4:v0 = p0 = 0 **/
  for(k = 0; k < Msize; k++){ v[k] = 0; p[k] = 0; }  

  /** step5: BiCGSTAB begins  **/
  for (i = 1; i <= MAXIT; i++) {
    if((i % nverbose) == 0) printf("Iteration %5d : ", i);

    /* 5.1 */
    for(rho_1 = 0.0, k = 0; k < Msize; k++) rho_1 += conj(r0[k])*r[k];

    /* 5.2 */
    beta = (rho_1/rho_0)*(alpha/w);
    if(cabs(rho_1)==0 || cabs(rho_0)==0 || cabs(w)==0){ 
      printf("\t Division by Zero in 5.2 \n"); 
      break;
    }
    
    /* 5.3  new p */
    for(k = 0; k < Msize; k++) p[k] = r[k];
    if (i > 1) {
      for(k = 0; k < Msize; k++) 
	p[k] += beta*(p[k] - w*v[k]);
    }

    /* 5.4 preconditioned phat and A*phat  */
    for(k = 0; k < Msize; k++) phat[k] = p[k];
    for(k = 0; k < Msize; k++) // (S - w^2) * p_hat
      v[k] = S1[k]*phat[k] + S0[k];
    //      v[k] = orb_erg[k]*phat[k] + potential[k] - k2*phat[k];
    
    /* 5.5 new residual and then alpha  */
    for(temp = 0.0, k = 0; k < Msize; k++) temp += conj(r0[k])*v[k];
    if(cabs(temp) == 0){
      printf("\t Divide by zero.\n");
      break;
    }
    alpha = rho_1/temp;

    /* 5.6 new direction s out of v  */
    for(temp = 0.0, k = 0; k < Msize; k++) {
      s[k] = r[k] - alpha * v[k];
      temp += conj(s[k])*s[k];
    }
    resid = sqrt(cabs(temp))/normb;
    if((i % nverbose) == 0) printf("residual = %15.10lf,", resid);
    if(resid < TOL) {
      for(k = 0; k < Msize; k++) x[k] += alpha * phat[k];
      printf(" and it is small enough to finish BiCGStab iteration.\n");
      free(r); free(r0); free(p); free(phat); 
      free(s); free(shat); free(t); free(v);
      return 0;
    }

    /* 5.7 precondition s to shat and calculating the t = A * shat */
    for(k = 0; k < Msize; k++) shat[k] = s[k];
    for(k = 0; k < Msize; k++) 
      t[k] = S1[k]*shat[k] + S0[k];
      //      t[k] = orb_erg[k]*shat[k] + potential[k] - k2*shat[k];

    /* 5.8  */
    temp1 = 0.0; temp2 = 0.0;
    for(k = 0; k < Msize; k++){
      temp1 += conj(t[k])*s[k];
      temp2 += conj(t[k])*t[k];
    }
    w = temp1/temp2;

    /* 5.9 obtain the new solution x */
    for(k = 0; k < Msize; k++) x[k] += alpha * phat[k] + w * shat[k];

    // 5.10 convergece test
    for(temp = 0.0, k = 0; k < Msize; k++){
      temp1 = alpha * phat[k] + w * shat[k];
      temp += conj(temp1)*temp1;
    }
    normr = sqrt(cabs(temp));
    resid = normr/normb;

    if((i % nverbose) == 0) printf(" delta_x = %15.10lf.\n", resid);
    if (resid < TOL){
      break;
    }

    /* 5.11 calculate new residual r */
    for(k = 0; k < Msize; k++)
      r[k] = s[k] - w * t[k];

    // and prepare for the next iteration
    rho_0 = rho_1;
  }

  free(r); free(r0); free(p); free(phat); 
  free(s); free(shat); free(t); free(v);

  return 1;

}
